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ABSTRACT 

The numerical solution of the full Navier-Stokes Equations for 
viscous flows with high Mach numbers and a strong detached bow shock 
is obtained. Two dimensional flows around (a) a circular cylinder, and 
(b) a circular cylinder with an aft-body in the form of a fairing, have 
been considered. The solution of the compressible M.S. Equations was 
accomplished by the method of finite differences. An implicit scheme 
of solution, the S.O.R. , was used with the optimim acceleration parame- 
ters determined by trial and error. The tensor notation was used in 
writing the N-S Equations transformed into general curvilinear coordi- 
nates. The coordinate system used is a general non-orthogonal 
curvilinear system with coordinate lines coincident with all boundaries. 
This coordinate system is numerically generated. The computational 
domain is limited upstream by a boundary located at a short distance 
ahead of the bow shock. T;.e shock has been treated as a sharp but 
continuous transition zone. (Shock capturing method.) The boundary 
conditions at this upstream boundary are the uniform flow conditions. 

The Euler equations were solved on the exit plane to establish the 
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downstream boundary conditions. The capability of attracting the 
coordinate lines to other p re-designated lines and, or, points exists. 
This technique was applied to blunt body flows with strong-shock to 
concentrate in and define the region of the shock. The amount of 
concentration was controlled by the factors "attraction amplitude" and 
"exponential decay" which were expressed as functions of the local 
density gradient across the shock. The equations for the generation of 
the coordinate system (using coordinate control) were solved, followed 
by the solution of the N. S . Equations , at the end of a set of given 
number of time steps. This process was repeated for different sets; 
thus the coordinate system concentrated in the shock region and moved 
with the shock. "Wiggles", by far constituted the one major problem 
that needed to be overcome. These oscillations that were encountered 
had no physical meaning and gave rise to quantities such as negative 
temperatures, which ultimately caused the computational program to 
break down. It was found that the application of certain dissipative 
finite-difference schemes damped these oscillations. The Shuman Filter, 
in particular, proved the most effective. The results are obtained for 
a Mach No. of 4.6 and a Reynolds No. of 10,000, and are compared with 
the experimental results available for these free stream conditions of 
flow around a circular cylinder. 
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CHAPTER 1 


INTRODUCTION 

In the accurate numerical prediction of real flows in general, 
there remain serious difficulties due to the complexity of these flows 
arising from unsteadiness, turbulence, three-dimensional character, 
large variations of flow gradients, etc. The numerical solution of the 
full equations for viscous flows, i.e., the Navier-Stokes Equations has, 
thus, attracted much attention, to treat fluid flow problems for which 
no simplified model exists. Generally it is found that the N-S 
Equations apply when the continuxim hypothesis is satisfied; this 
condition being satisfied when the Knudsen Number is less than 0.01. 
Victoria and Widhopf [1] reported that the use of the N-S Equations to 
solve hypersonic low Reynolds/high Knudsen Number flow problems, is 
valid. They considered the flow about a sphere in a Mach 10 flow at a 
free stream Reynolds Number of 152, with a corresponding Knudsen Number 
of 0.10. The range of validity of the N-S Equations does, indeed, seem 
to cover most of the aerodynamic problems relating to aeronautics and 
astronautics. 

The present investigation concerns the numerical solution of the 
complete N-S Equations for viscous flows, with strong shocks, ahead of 

a two-dimensional arbitrary body. The detached bow shock separates 
the flow disturbed by the body from the undistrubed flow. Behind the 
shock wave there is a sub-sonic region bounded by it, the body and the 
"sonic lines." The flow field outside of this sub-sonic region is again 
super-sonic, with the exception of the region of the wake immediately 
behind the body. The shape of the bow shock ahead of the body is 



Influenced by the shape of its leading edge. The numerical solution of 
the compressible N-S Equations in the case of a mixed sub-sonic and 
super-sonic flow field as described above is a difficult task and cannot 
yet be ranked along with the standard theoretical tools currently used 
in applied aerodynamics. 

Of the various types of methods in izse for the numerical solution 
of the N-S Equations, the finite-differencs method is perhaps, by far the 
most widely used. In the last few years, finite-element methods have 
received increasing attention as a substitute to finite-difference 
methods in fluid mechanics problems, particularly in the solution of 
the incompressible N-S Equations. 

To keep the computing time within reasonable bounds it is important 
to minimize the number of mesh points and this usually requires that 
the mesh system be taken non-uniform in the physical plane. This can 
be achieved by imposing a variable mesh spacing in a given coordinate 
system. The mesh system is generally also chosen so as to make the 
boundaries of the computational domain, particularly solid walls, 
coincide with mesh lines; this considerably simplifies the treatment 
of boundary conditions. One of the novel features of the present 
research is the use of the body fitted coordinates which are a system 
of general curvilinear coordinates, that are niimerically generated, 
developed by Thompson, et. al. , [ 2 ] and Thames [ 3 ]. This system of 
coordinate generation is quite versatile due to the fact that it allows 
the coordinate lines to be coincident with all boundaries of a general 
multiply connected region including the boundaries formed by the solid 
walls of any number of quite arbitrarily shaped bodies. The boundaries 
may even be time-dependent. Thus with this procedure the numerical 
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solution of the N-S Equations tnay be obtained on a fixed rectangular 
field or the transformed plane , with a square mesh. No interpolation 
of the flow variables is required regardless of the shape of the 
physical boundaries or the spacing of the curvilinear coordinate lines 
in the physical plane. 

Once the coordinate system has been chosen , and the choice of a 
mesh size is made, various techniques may be made use of to reduce the 
computing time as much as possible for a given numerical scheme. In 
an explicit scheme of solution, the local maximum time step depends 
strongly on the local mesh size in physical space; if the physical mesh 
varies considerably throughout the computational domain, the time step 
will be determined by the smallest mesh and will be, in all probability, 
very small. It is then practically indispensable to divide the domain 
into several regions in each of which a different time step is used so 
as to reduce the total number of operations necessary to advance the 
solution in time in the entire field - all of which leads to compu- 
tational complications and restrictions. Thus it is seen that time step 
limitation is the main drawback of explicit schemes; however, these 
methods have been widely used because of their simplicity and the fact 
that the number of numerical operations at each step is kept to a 
minimum. Another approach which is attracting more attention now is 


the use of implicit schemes which lead to less severe stability con- 
ditions or which are unconditionally stable. Roger Peyret and 
Henri Viviand [ 4 ] concluded (in 1975) that *’no clear cut conclusion 
can be drawn at this time regarding the best type of method - implicit 


or explicit." However, the fact that an implicit scheme is uncon- 
ditionally stable proved irresistable and one such scheme - the S.O.R. 
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was used in the present work. Also tensor notation was used in writing 
the N-S Equations transformed into general curvilinear coordinates as 
this measure greatly simplified the overall problem of analytical 
development and rendered the method more general. Thus the extension 
of the method to three dimensions is now purely formal. 

The flow considered being super-sonic about a blunt body, a 
detached bow shock exists ahead of the body, as noted earlier. The 
computational domain was limited upstream by a boundary located at a 
short distance ahead of the bow shock, treating the shock as a thin 
but continuous transition zone. This is the "shock capturing" method. 
Tannehill and Holst [ 5 ] used this approach for low Reynolds Number 
flows. However, they found that (as indeed, we too) for high Reynolds 
Number flows it was not "practical" to capture the bow shock because of 
the numerical difficulties associated with the large gradients at the 
bow shock. Instead, they foxmd it more convenient to treat the bow 
shock as a discontinuity, across which the Rankine-Hugoniot Equations 
could be applied, while leaving the boundary layer to be "captured." 
This approach is called the "shock-fitting" method and is generally 
favored because, the flow being practically inviscid in the vicinity 
of the shock, the N-S Equations are not really needed to calculate the 
shock; but the problems associated with this method are those that arise 
due to the necessity of having to couple the inviscid and viscous flow 
solutions and the a priori definition of inviscid and viscous regions. 

The shock-capturing method was used in this research, not withstanding 
the numerical difficulties associated with this method, because it is 
much more convenient to solve the N-S Equations in the entire flow field 
when the inviscid flow region is of a small extent. This situation 
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exists in the case of bow shocks in high Mach number flows. Now the 

thickness 6 of a normal shock wave in real viscous fluids varies as 
s 

1/Ee for fixed Pr. For high Reynolds Number flows it then becomes 

likely that 6 < Ax (the mesh size normal to the shock). When this 

s 

occursj oscillations develop aft of the shock leading to the numerical 
difficulties as reported by Tannehill and Holst and a host of other 
investigators. If 1/Re 0 and if the finite-difference method has no 
artificial viscosity, then there is no dissipative mechanism by which 
the oscillations may be damped out. Many authors have associated these 
spatial oscillations, also called "wiggles," with non-linearities, 



or with linear instabilities in the transient calculation. Roache [ 6 ] 
has however, demonstrated that wiggles are not caused by iterative 
instability or non-linearities, but that they simply are a solution 
of the finite-difference equation used. In the ntunerical computation 
of flow fields containing shock waves, these oscillations or wiggles 
may be damped by the application of dissipative finite-difference 
schemes. The concept of introducing an implicit dissipation term by 
using a dissipative finite-difference scheme to damp out short wave 
oscillations has been used by many researchers. Vliegenthart [ 7 ] 
reported that these encountered oscillations which have no physical 
meaning, can be suppressed by applying Shuman's technique of introducing 
dissipation, and that, in certain cases this even appears to remove 
nonlinear instabilities as in the case of the computation of a detached 
shock in front of a blunt-body. Due to its extreme simplicity and 
effectiveness the Shuman Filter, which is discussed in Chapter IV, 
has been used in the present research, to overcome the problem of 
wiggles. Unlike other schemes, the Lax-Wendroff, for instance, the 
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Shumann Filter is not an integral part of the difference approximation. 
Lax's scheme, yields a shock which is spread out over a large number 
of mesh spacings, while the Shumann Filter yields a shock which is 
considerably narrower. 

The capability of attracting the coordinate lines to other pre- 
designated coordinate lines or grid points exists at present, and 
through the application of this technique to blunt body flows with 
strong shocks, an effort was made to concentrate in, and define, the 
region of the shock. The magnitude of concentration is controlled by 
the factors "attraction amplitude" and "exponential decay." (These 
terms are discussed in Chapter II.) These cooordinate control factors 
were expressed as functions of the local density gradient across the 
shock. The equations for the generation of the coordinate system, 
using coordinate control, were solved, followed by the solution of the 
N-S Equations, at the end of a sat of given number of time steps. This 
process was repeated for different sets of time steps. Thus the 
coordinate system was refined in the region of the shock and moved with 
the shock. This measure obviated the need of having a very refined 
mesh in the entire computational domain by providing refinement only in 
the region through which the shock happened to be passing at any given 
time. 

The results presented pertain to a flow about a two-dimensional 
circular cylinder with a free stream Mach Number of 4.6 and a Reynolds 
Number of 10*^. A few results for the case of a circular cylinder with 
an aft body (fairing) have also been Included for comparisons. The 
wall pressure, normalized with respect to the stagnation pressure, as 
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reported by TannehlU and Holst [ 5 ] has been compared with the present 
solution for the circular cylinder and the two results are in good 
agreement. All the numerical results are discussed in Chapter v. 
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CHAPTER II 


NUMERICAL GENERATION OF CURVILINEAR COORDINATES 

2.1. The Botmdarv-Fltted Coordinate System 

The accurate numerical representation of boundary conditions, be 
it the body surface or the infinity boundary, is best accomplished 
when the boundary is such that it is coincident with some coordinate 
line. Under such circumstances the boiandary may be made to pass 
through the points of a finite difference grid constructed on the 
coordinate lines. This eliminates the need for any interpolation of 
the dependent variables between points of the grid. The avoiding of 
interpolation is particularly important for boundaries with strong 
curvature or slope discontinuities, and for differential systems that 
produce large gradients in the vicinity of the boundaries. Thus the 
generation of a curvilinear coordinate system with coordinate lines 
coincident with all boundaries is an essential part of a numerical 
solution. Extensive use of the method of numerical generation of a 
curvilinear coordinate system due to Thompson, Thames and Mastin [ 2 ] , 
has been made in the present work. The main idea of this method is to 
fill the computational domain enclosed between the body and the 
external boundary with Intersecting coordinate lines in the physical 
(x,y) space. 

Let C = ?(x,y) and n = n(x,y) be two continuously differentiable 
functions of the Cartesian coordinates (x,y). Further, let n = = 

constant be the body contour, while n “ “ constant be the external 

boundary contour. The region S n ^ must now be filled by inter- 

secting coordinate curves € = constant and n = constant. Because of 


^1 


^1 








the closed region under consideration it is natural to specify the 


determining differential equations for ^ and n as elliptic equations 


to he solved under the proper boundary conditions for | and n at the 


body and at the external boundary. Since the simplest elliptic equation 


is the Laplace Equation, we then pose the problem of solving the 


Laplace Equations for 5 and q with x and y as independent variables 


under the Dirichlet boundary conditions. Let be the curve 


defining the body contour q “ q and F„ be the curve defining the outer 

o ^ 


boundary q => q in the xy-plane as shown in Flg.l. The elliptic 


boundary value problem is then 


7^5 = 0 


( 2 . 1 ) 


V^q = 0 


( 2 . 2 ) 


on Fj^: C = f^(x,y), n = 


(2.3) 


on F 2 : C = f„(x,y), q = q^ 


(2.4) 


The solutions of Eqs. (2.1) and (2.2) under the boundary conditions 


(2.3) and (2.4) can conveniently be obtained in those cases when q^ and 


q can be specified by simple analytic methods (such as a circle. 


ellipse, etc.). To obtain coordinates for arbitrary shaped bodies, it 


is convenient to transform the Eqs. (2,1) and (2.2) such that x and y 


are the dependent while 4 and q are the Independent variables. This 


transformation is more easily performed for either two or three- 


dimensional coordinates by the method of tensor analysis and is detailed 


in Appendix B. Referring to Eqs. (B-13) and (B-14), we find that 


Eqs. (2.1) and (2.2) are equivalent to 


hi hi ■ ^hi hn * “ 


C2.5) 


hi hi - ®11 


( 2 . 6 ) 






where the variable subscripts denote partial differentiations and the 
is the metric tensor defined in Appendices A and B. The boundary 


conditions are now 



On rj: 

X » F^(5,n^), y « F2(5,n^) 

(2.7) 

On F*: 

X “ y - 62(5,nJ 

(2.8) 


where as shown in Fig. 2, F* and F^ are the images of the body and the 
external boundary contours in the 5n“plane. 

The geometrical meaning of the transformed equations (2.5) and 
(2.6) is that the body and the external boundary contours in the xy- 
plane have been mapped on to the Cn“Plane which is rectangular. In 
other words, we can say, that the contours in the xy-plane have been 
opened up to form the straight lines n “ >1^ * constant, and n “ " 

constant in the 5n-plane. This can be achieved by imagining a cut 
connecting the body and the external boundary in the xy-plane as shown 
in Fig. 1., such that all functions and their derivatives are con- 
tinuous in crossing the cut. Since a cut line is a part of the field, 
no boundary conditions can be imposed on F^ and F^ of Fig. 2. 

The appearance of n and n in Eqs. (2.7) and (2.8) is now purely 

o ^ 

symbolic, denoting the names of the body and of the external boundary 
respectively. Given the body and the external boundary contours, we 
can always establish the values of x and y either graphically or 
analytically for any desired distribution of 5-values. The n-values 
can be chosen arbitrarily to form rectangular meshes in the 5n-plane. 

Equations (2.5) and (2.6) are the basic equations for the 
generation of coordinates. To have a control over the spacing of the 
5 and n lines, xje envisage another general transformation, say from 
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CjH to ?■* and n''. Retaining the 5>n notation, the equations take the 
form 


822 ■ ^^12 =^511 «u *tm ■ **5 

S 22 ^{5 - 2Si 2 ^Cn * *11 ■*• % 


For details on the above derivations refer to [ 8 ] . The same form of 
equations are obtained if one starts considering the Poisson Equations 
in place of Eqs. (2.1) and (2.2) and inverts the transformation from 
x,y to C»n as independent variables. 

The function P and Q are to some extent arbitrary and can be chosen 

in various ways to have a desired distribution of coordinates in a given 

region. In the present research we have made P and Q to depend on the 

density gradients to contract the coordinates in the region of the 

shock. The chosen forms of P,Q are [ 9 ] 

_ n' 

P - g ^2^ Aj; sgn (C-5j^) exp 


m 


+ g 2 B' sgn (5~?^) exp (-E£Rj^) 


Q = g 2 sgn (n-tlj^) exp (-Dj^ln-Hjj^l) 
h=l 


( 2 . 11 ) 


+ g 2_ sgn (n-Hjj) exp (-Ej^R^) (2,12) 

i— 1 

where 

g “ 811 822 “ <812^^ ’ 

h “ (2.14) 

The first terms on the right hand sides of both (2.11) and 2.12) are 
used in the line attraction, while the second terms in both equations 
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are used for the point-attraction. Various terms which appear in these 
equations have been defined in the "List of Symbols". 

In the present research, we have used only the point attraction 
term of Eq. (2.12) to concentrate the coordinate lines near the shock. 
The amplitude factor is a function of the maximum density difference 
along each C-line in the region of the shock. This amplitude factor 
thus changes according to the position (C,n) and is defined as 

= (constant) (p2“Px^^‘^l (2.15) 

where the subscripts 1 and 2 denote the respective values in the front 
and behind of that shock which has been computed without coordinate 
contraction. The constant appearing in (2.15) is selected by trial 
and error and retains the same value for all ^ and n positions, that 
correspond to the shock location. 

The method of numerical coordinate generation offers much freedom 
in the orientation of both the ^ and n coordinates in the physical xy- 
plane. For example, the n = const, lines can be chosen to go round 
the body as shown in Fig. 1, or they may not be chosen to form a 
complete circuit as shown in Fig, 5. However, a suitable choice has 
to be made in advance of computing the coordine *:es , because the 
resulting configuration of the body segment, the cut lines, or the re- 
entrant segments, and the outer boundary segments in the transformed 
5n-plane depend on this choice. In the present research we have chosen 
the coordinate configuration as shown in Fig. 3, in which the front 
outer boundary is a hyperbolic arc and the rear outer boundary is a 
circular arc. Figure 6 shows the corresponding segment orientation in 
the plane. This t 3 fpe of segments orientation requires much care in 
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the computer programing for both the coordinate generation and for the 
numerical solution of the Navier-Stokes Equations. In Section 2.2 the 
finite-difference approximations of all the equations have been 
discussed. 

2.2. Finite Difference Approximation 

In this section the finite difference approximations and the 
numerical methods used for the solution of the elliptic systems, Eqs. 
(2.9) - (2.10) is discussed. 

Before proceeding mith the pertinent method of solution, it is 
important to mention that in the method of body-fitted coordinates it 
is superfluous to specify the step sizes AC and Ap, both of which are 
equal to 1. If DIAX and Jl'lAX represent integers for the maximum 
numbers of C and rj points respectively in a field, then this input and 
the desired contraction controlled by the amplitude and decay factors 
of Eqs. (2.11) - (2.12) decide the variable mesh sizes to be obtained 
by solving the generating system (2.9) - (2.10). This aspect has been 
thoroughly discussed in [ 8 ] . Thus the main utility of ntimerically 
generated body-fitted coordinates actually lies in the availability of 
meshes or nets in the Cn-plane on which the Navier-Stokes Equations are 
to be solved without specifying the step sizes. Further, the vari- 
ations both along the C-and n-coordlnates , are labeled by the con- 
secutive Integers in the range 1 I < IHAX and 1 < J S JMAX. 

The solutions of Eqs. (2.9) - (2.10) have been obtained by the 
Gauss-Sledel method with successive over relaxation (SOR) under the 
prescribed boundary values for x and y on the body and the external 
boundary contours, along with the prescribed values of UlAX and Jl'lAX. 


13 


The spatial derivatives are approximated by the central differences 




(2.16) 




(2.17) 


where X is a surrogate variable and I»J denote positions along ^ and n 
coordinates respectively. Similar expressions are obtained for X , X 

n nn 

and X . 

5T) 

The solutions of Eqs, (2.9) - (2.10) yield x and y for the whole 

flow field as functions of ^ and q. This data is then used to generate 

the derivatives x^, x , y^, y * the metrical coefficients g, the 
C’ n’ ■’n ij 

i 

determinant g, and the Chrlstoffel symbols 

It was mentioned in Section 2.1. that the orientation of the 
coordinates of the form shown in Fig. 5 requires due care in obtaining 
derivatives on the cut line. Referring to the schematic shoxsn in Fig. 

3 with the x-axis oriented along the cut line, let 1^^ and denote 
the integral values of I on the upper and lower parts of the cut 
respectively. Thus 




(2.18) 


Equation (2.18) establishes the following correspondence between I. 


and lyj.; 


1 Corresponds to 


IMAX 


IMAX-1 


IMAX-I+1 


(2.19) 
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where and represent the same point of the body reached by t e 
lower and upper parts of the cut respectively. Obviously 

= (IMS + 1) - Igi . 


From (2.19), we conclude that 




( 2 . 20 ) 


x(I^^,+l,J) = x(Ij^jj-l,J) ^2.21) 

The first derivatives on the lower and upper parts of the cut are 




^yn^Lc " I 


where 2 < < Igi “ 1 • 


(VUC ” 2 

cy^^uc “ i " y(iuc-i,i)i 




( 2 . 22 ) 


(2.23) 


where 


I„ + 1 £ I„C S ^ 








CHAPTER III 


THE NAVIER-STOKES EQUATIONS 


3.1. Formulation of the Problem 

For solving the blunt body problem in general curvilinear coordi- 


nate system, we consider the nondimensionalized Navier-Stokes system of 
equations in the invariant tensor form. The conservation equations are 
Nass Conservation: 


^ + div(py) - 0 


(3.1) 


Momentum Conservation: 


(py) + div T = 0 


(3.2) 


Energy Conservation: 


+ div b = 0 


(3.3) 


where 


T ~ pyy + pi - so 


b =■ (’F+p)y - ea»y - 


OqP grad T 


f - pa + I- p|yp 


0 « K1 + d 


K = A div V 
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d » wdef y - w [grad y + (grad y) ] 

e « 1/R 
e 


“O ■ 


<3.4: 


The nondiinensionallzed equations (3.1) - (3.4) have been obtained 
by referring all lengths to the diameter 2R*j velocity vector, density, 
viscosity, temperature, enthalpy and pressure to the free stream values 
V*, p*, p*, T*, h*, and p* V*^ respectively. The N-S Equations (3.1) - 
(3.4) by themselves do not give a complete description of the motion of 
a compressible fluid because changes in pressure and density cause 
temperature variations and thermodynamic principles must, therefore, 
enter into the considerations. Asstmilng the gas to be caloricaily and 
thermally perfect, the equation of state in the nondimensional 
variables is 


P ■ 

where y is the adiabatic exponent, 
temperature is given by 


pL 

Similarly, the nondimensional 


(3.5] 


X « y(y-1)m 2 (1 _ I |yl2) 


(3.6) 


The relation between temperature and viscosity is provided by the 
Sutherland formula, which in nondimensional form is 



V - 


(l+Sjj)T^^^ 

-T+-sf 


s! = 110“K . 


(3.7) 


\^heve S, 





Since ¥ is a function of T, hence the system of equations (3.1) - (3,7) 
form a closed system of equations provided X = X(p) is also prescribed. 
In the present formulation we have used the Stokes' condition 

3A + 2y = 0 (3.8) 


as the required relation. 

The boundary conditions for the system of equations (3.1) - (3.3) 

t t 3T 

at the body surface are: jy| « 0, T = T^, or (g|^)^ specified 


P T 

- w w 
Pw" 


(3,9) 


at free stream infinity; |v[ 


P = 


1, p = 1, T = 1 
1 


¥ =j + 


Y(y-1)m2 


In (3,9) the subscript w denotes the wall condition. The density is 
not known in advance but must be obtained by the equations themselves. 

Since the governing equations are of parabolic-elliptic type, we 
therefore need to specify the outflow boundary conditions. In place of 
specifying the derivative conditions, we have used the complete 
solution of the Euler Equations to specify the downstream boundary 
condition with the stipulation that at the outer downstream boundary 
the effect of viscosity is negligible. 
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3.2, Transformation of the Equations 

The governing equations (3.1) - (3.3) are now transformed to a 

i 

general coordinate system 5 (i = 1»2) by the method of tensor analysis 
(Appendix A). Using the summation convention on repeated indices, and 
introducing the new dependent variables 

o " p, P “ t^g” p, , E « (3.10) 

the equations are 

^ + -4r (ov"'-) » 0 (3.11) 


(ovS +-\ (avV) “}■ ^ 


3r 


■jk 


ar 


- P r|. g^^ + ei-K iK'/r 


ag' 


(3.12) 


+ {(E + P)v^} * e — ^ (*^ v^) + Of. (w*^ “^) (3.13) 

f,..* ....a* u ..i-i acJ 


3E L a 
at 


ag* ag* '' ag 

where ®ij is the transformation metric tensor, and 




ag** 
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(3.14) 




■Vj -.■J.- 

‘>1 

j ‘ i ' i : '■ 


i L j V V;.- 




-.v'l ' ■■! -rf ■•'.I 


-.■■■!• 



In Eqs. (3.10) - (3.14), the superscript indices denote the 
contravariant components, while a comma denotes the covariant 
differentiation, i.e. , 


sjJ 3fc 


where rt, is the Christoff el symbol, 
jfc 


(3.15) 


Equations (3.11) - (3.13) are applicable both to two and three 
space dimensions. In two dimensions, writing for brevity 
5^ ** « n, v^ “ u, *» V 

A « ou, B “ ov, H “ e(u)^, M = e(v)^» N » ouv, (3.16) 

Q » (E+P)u, R « (E+P)v 
the equations become of the form 


la + M + M „ 0 

3t 35 3n 


(3.17) 


f f I? + ’’iiS + ■'22» - ® 


(3.18) 


~ + — ■}• “ + r® H + 2F^ N + r^nM * i 
3t 35 3n 11 12 22 ^ 


(3.19) 


3t 35 3n ^ 


( 3 . 20 ) 


The expanded form of various terms appearing in Equations (3.17) - 
(3.20) are given in Appendix A. 

The boundary conditions for Equations (3.17) *- (3.20) are at the 

3T 

body surface: u = v = 0, T = T^, or specified 


E 


w 


T a , 
w w 

Y(y-1)m2 


\ “ (Y-l)E^ . 
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at free stream infinity: u » v ~y^//g 

0 ^ 

T "»• 1 


P -> 1^/ymJ 
E *»■ + 
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(3.22) 


where a ^ or n subscript implies partial differentiation. 

The local Mach Number Mjj^ is given by 

® (3.23] 

where 

lyl^ * (3.24] 

and 

T . 'rh-DH* (f - I lyl*) (3-25: 

The relations between the local Cartesian and the local contravariant 

components of the velocity vector y are 

U » ux„ + vx 

5 n 

V “ ’ly^ + vy^ . (3.26] 


Tlie vorticity o) is given by the formula 


to 


i ,!ai 


3v, 






(3.27] 


where v^^ are the covariant components of y, which are related with the 


contravariant components as 
v.< 




Su “ + «12 » 


(3.28] 


The pressure coefficient is defined as 


p‘-p; 


i. n* v*2 


CO 


c 

p 


(3.29 



CHAPTER IV 




i. ij- 
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We now discretize Eq. ( 4.1 ) by a fully implicit difference 
approximation. The time derivative is approximated by a first-order 
backward difference at n-M, where n is the time step of size At, while 
the spatial derivatives are approximated by central differences. The 
system of equations are solved by polnt-SOR, which are 

(4.6) 


n+1 n+lCp) , _ 

Wt t ” Wt t + 


::i,J "I,J • '“i'l.J 
where tu is the relaxation parameter and the superscript (p) denotes 
values at the previous iteration* The fxmction R is 


~I,J 


n n-i-lCp) ^ TjU+1 

I,J ” ~I,J 2 ^-I+1,J “ “I-1,J 


. pU+1 _n+l . 




(4.7) 


where the values of w in F, G and H are those which are the most recent 
values available from the previous Iteration. Fully expanded forms of 
(4*7) are given in Appendix D. 


4.2, n-Derivatives on the Cut 

To find the q-derivatives on the cut, we refer to Figure 4. The 
point 1 of the physical plane (Figure 3) transforms to the location 
marked x on the transformed plane, while point 2 remains at its original 
position relative to the cut. Thus in principle, the function value at 
the fictitious point shown under the lower cut must be replaced by 
the function value at the point 2 on the upper cut. Now two cases 
arise depending on whether the function is a scalar or a vector. 
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Tj“derivatives are 


For a scalar f, the first 


2 

- 2 - *<ioc>2>l 


/•3f. 1 


Thus 


(#). 


-(#). 


Similarly, the second n-derivatives are 


(4.8) 


(^) 


3n2 ' 2f(Ij.^,l) + 


(4.9) 


(^). 


(4.10) 


arj2 “ 2f(I^^,l) + f(Ij^^,2) 

But f is a scalar, so that 

hence both (4,9) and (4.10) represent the same value. 

lo find the n-derlvatlves of a directional quantity u on the cut, 

we need the value of u at the fictitious point. Since In the physical 

Plhne u in the lover part of the cut la directed opposite to that on the 
upper part of the cut, hence 


9Uv 


“2 f"(^c*2) - (-u(Iy^,2))] 


me (iU) 


(4.11) 


The same holds for v, and 

Baaed on the preceding analysis It is easy to shoe that either for 

a scalar function f or a vector function y the derivatives across the 
cut U 2 Te contXuuouSa 



The wall density is calculated by evaluating each term of the 
continuity equation (3.17) at the body surface. Denoting by J=1 the 


body surface, we have 




Using a three-point fonrard difference approximation for the right hand 
side, we obtain 

(4.12) 


°I,1 “ “"l,! " 2 “ ®I,3^ 


where = 0. 

igX 


Though Eq. (4*12) is fully implicit, nevertheless its use at the 
trailing edge point always produces unrealistic density values. To 
circumvent this difficulty, Eq. (4.12) was used at all points of the 
body except at the trailing edge point where an explicit scheme based 
on the leap-frog method was used: 

n 


n 


(i£)“ , . (3B) 


n+1 _ n-1 

f.l ~ I»1 _ 1 /«n . 

2At 2 ^®I,2 " ®I,0^ 


where ”0" is a fictitious point. Using a second order extrapolation in 
space and time, we get 


.n 


n" _ 2 r®""^ _ R^~^ 

^I.O ®I,2 


Thus the wall density is obtained by the expression 
,n 1 1 .n A ±- /'r.n i ^^n— 2^ 


„nrl n ... /^n . 

Ot 7 = o-r 1 - 4t (B.^ ., + ,) 


(4.13) 




4.4. Nonlinear Instability 

In the solution of compressible flow equations, several types of 
nonlinear instabilities are encountered. Among these, the most dominant 
is due to the difference approximation of the convective derivative. 
These instabilities can be avoided by introducing some dissipation in 
the difference approximation of the differential equation being solved. 
In this context, McCormack [10] has used a fourth^order damping term 
for his explicit schemes. Boris and Brook [11] have developed a flux- 
corrected transport (FCT) technique idiich is quite efficient for use 
in the continuity equation. Vliegenthart [7] and Harten and Zwas [12] 
have iised "Shuman Filtering" to supress the convective instabilities. 

In the present research both the FCT of Ref. [11] and the Shuman 
Filtering of Ref. [7] have been tried. Though^ the Shuman Filter adds 
more dissipation than desired, particularly near the shock, it always 
produced wiggle-free solutions for all regions of the flow field. The 
application of Shuman Filter amounts to replacing the flow variable 
w? , by w? . = w^ in 4.7, where 


i,j 




(4.14) 


This scheme was used on all the four primary flow variables in the 
form of a/V^, kf/gt B//g and E/v^. After Equation (4.14) was applied 
to these quantities the original dependent variables of the N-S 
Equations taere recovered by multiplying each of the filtered variables 
by the Jacobian, This filtering technique was carried out on the 

converged solution that was obtained at the end of 5 time steps, each 
time step (At) being 0.01. The frequency of application of the Shuman 
Filter had to be determined by trial and error. 
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4.5. Downstream Boundary Conditions 

For the solution of the parabolic-elliptic system of Equations 
(3,17) - (3.20), beside the boundary conditions at the body surface and 
at the upstream free boundary (Eqs. 3.22)), it is also necessary to 
specify a proper set of conditions at the doxmstream boundary. To 
obtain these conditions, the downstream boundary was placed at a 
sufficient distance where the viscous dissipation is negligibly small 
and the complete Euler’s Equations were solved on the entire down- 
stream boundary. The computer program has been structured in such a 
way that it solves both the Navier-Stokes and Euler's Equations with 
each iteration. 

The Euler Equations that were solved on the downstream boundary, 
at each iteration, are: 

(4.15) 


80 , 3A . SB _ 




(l/g)[gi2 an ~ §22 3£^ ^12^ 


’12 3n ®22 3£ 


“ ®12 ^^12 ^ 22 ^^ 


(4.16) 


If If H ^ 


11 


12 ‘ 


22 


- (l/g)[gj^2 3£ ■* ®11 9n^ ^ (p/g)tSl2.^^12 ^22^ 


" %2 ^^11 *** ^ 12 ^^ 


9E + ia+ 3R 0 

3t 3£ 3n 


(4.17) 

(4.18) 
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The finite-differencing of these equations was carried out by 

three-point forward-differencing in the £ direction on the C = 1 

coordinate and by three-point backward-differencing in the C direction 

on the 5 = C coordinate. Central differences were used for the n 
max 

derivatives. 

^•6. Coordinate Contraction Hear a Shock 

The capability of attracting the coordinate lines to other pre- 
designated coordinate lines or grid points exists at present and 
through the application of this technique to blunt body flows with 
strong shocks an effort was made to concentrate in, and define, the 
region of the shock. The magnitude of concentration is controlled 
by the factors B. and defined in Eq. (2.12). These coordinate 
control factors were expressed as functions of the local density 
gradients across the shock as shown in Eq. (2.15). The equations for 
the generation of the coordinate system (Eqs. (2.9), (2.10), using 
coordinate control, were solved, as well as the Navier-Stokes 
Equations, when a quasi steady-state has been reached. This process 
was repeated after a pre-assigned number of time steps. Thus the 
coordinate system was refined in the region of the shock and moved 
with it. This measure reduced the need of having a very refined mesh 
in the entire computational domain by providing refinement only in 
the region through which the shock happened to be passing at any 
given time. However, it must be noted that this refinement near the 
shock was achieved at the expense of the accuracy near the wall where 
a finer mesh is always needed to resolve the boundary layer. 
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CHAPTER V 

COMPUTATIONAL PROCEDURE AND RESULTS 
5 .I 4 Computational Procedure 

The first step in the ntunerical solution of the transfoinned N-S 
Equations Is the determination of the computational domain so that the 
appropriate boundary conditions may be prescribed arotind it. In the ' 
present case of super-sonic flow, the flow field remains unperturbed 
upstream of the bow shock wave^ and the computational domain Is limited 
upstream by a boundary located at a short distance ahead of the bow 
shock. The **stand-off” distance of this detached shock, for a given 
free stream Mach Number, is estimated using empirical equations [13]. 

On the upstream botmdary the uniform flow conditions are used as 
boundary conditions. Particular care had to be taken to ensure that 
the bow shock did not cut across any segment of this upstream boundary. 
On the downstream boundary the boundary conditions varied \rf.th time 
and were determined by solving the Euler’s Equations (cf. Chapter IV). 

The computational domain and the profile of the body in it having 
been determined, the next step was the numerical generation of the 
coordinate system which has already been described. The cartesian 
coordinates of each of the mesh points in the entire computational 
domain having been determined and stored, the coefficients that occur 
in the Navier-Stokes Equaitons due to transforming them into general 
curvilinear coordinates could now be calculated and stored in a file. 

The actual solution of the N-S Equations now starts with an 
assumed Initial guess of the solution for the entire computational 
domain. These initial conditions need not necessarily be physically 
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realistic and when they are notj the transient solution has no physical 
meaning. In the present case the initial conditions chosen for the 
whole flow field were the uniform flow conditions that were prescribed 
on the upstream boundary. It was however found that this could not be 
done if the free stream Mach Number was very high, or if the isothermal 
temperature condition prescribed on the body was far different from the 
free stream value. The finite- difference scheme chosen was the S.O.R. 
which is an implicit scheme. The value of the optimum acceleration 
parameter for all the equations, l.e., continuity, momentum and energy, 
was determined, by trial, to be 0.9. In approximating the convective 
derivatives of these equations, the average of the product (A.O.P.) 
finite-difference scheme rather than the product of the average (P.O.A.) 
proved fruitful, even though it is generally considered that a non-linear 
instability can result in regions of flow reversal when the average of 
the product scheme is used. 

The problem of the treatment of boundary conditions at an imperme- 
able wall in viscous compressible flows reduces to that of the 
calculation of the pressure or of the density. In this research the 
wall density was calculated from the continuity equation written at 
the wall. Peyret and Viviand [4 ] report that such a technique is of 
delicate use and may lead to strong oscillations or even to divergence 
if no artificial viscosity term is added to the continuity equation; 
and that, in particular, in the case of separated flows negative values 
of the density may be obtained. This, in fact, was what happened at 
the trailing-edge poinc using the continuity equation. This problem 
was overcome by using an explicit discretization based on the leap-frog 
scheme, only at the trailing edge (cf. Eq. 4.13). 
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la general it was found that at a given time-step iterative con- 
vergence to the tolerance of 10*“® occurred in about 6 iterations. 

T^hile carrying out these iterations the downstream boundary conditions 
varied with every iteration. Progress in time was made by increasing 
the time by a At of .01 at the end of each time step. The problem 
of wiggles was overcome by applying the Shumann filter at the end of 
every 5 time steps. Provision was made in the con^uter program to 
store the solution obtained at the end of any desired time step in 
a file; and the ability to read back from this storage file and to 
restart the program where it last left off was also incorporated. 

The coo^uter program also locates and calculates the maxiimim change 
that occurs in a typical flow variable, such as density, along every 
€ constant line that passes through the region of the shock. Thus 
at the end of any pre-designated time step the location of the shock 
and the change, in the density across it is automatically recorded. 

This information is used again if necessary in the generation of a 


new coordinate system wherein the coordinate attraction technique is 



used to refine the mesh in the immediate vicinity of the shock. 

On an average it took 0.525 minutes of computer time (on a UNIVAC 
1108) to achieve iterative convergence at each time step. The stand- 
off distance of the bow shock became quite constant after about 320 
time steps, each increment in the characteristic time-step being equal 
to .01. All the cases were, however, run up to 400 time steps 
the total computer time requirement to achieve this "steady-state" 
solution was about 2 hours and 30 minutes. 
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A reduction in the time-interval At should normally lead to more 
accurate results. The use of a filter, such as the Shuman filter, 
however, increases the dissipation with decreasing time-interval At. 
For the case of the clyinder, with a time-step of 0.005 the pressure 
distribution was generally found to be higher, particularly in the 
vicinity of the 90“ point, i.e,, the top of the cylinder. The use of 
the Flux Corrected Transport Filter was not quite satisfactory as it 
Introduced too little dissipation, as opposed to the Shuman Filter 
which introduced too much. In fact, the Shuman Filter introduced so 
much artificial viscosity that it overshadowed the effect of a 
reduction in the free stream Reynolds Number. For the same At of 
0.005 a reduction of Re from 100,000 to 10,000 made no difference to 
the solution. In general it was found that an increase in Reynolds 
Number required a smaller time-increment At. 



5.2. Discussion of Results 

The numerical solution of the complete Navier-Stolces Equations for 
a super-sonic flow was obtained for the flow about a two-dimensional 
circular cylinder. The uniform flow conditions used in the computations 
were M =* 4.6, Re = 10, 000 f T* = 167 °K, p* =* 14.93 N/M^ and P = 0.72. 
For these free stream conditions the coefficient of viscosity works out 
to be p* = 1.13154 x 10“^ kg/m-sec and the density p* = 3.11593 x lO""* 
kg/ra^. The ratio of the specific heats was assumed to be y = 1.40. 

All calculations presented in this report have been perfomed for the 
Isothermal wall temperature T* of 556°K. The diameter of the circular 
cylinder was 2R* = 0.3048m. The Knudsen Number for a perfect gas is 
defined by the expression i/yn/2 (M^/Re^), which for the above free 
stream conditions is 6.821529 x 10~**. 

The graphical results presented correspond to the steady state 
solutions at a characteristic time of 3.20, with the exception of the 
Mach contour plots which are presented at two different periods of 
time so as to give an insight into the formation and progress of the 
bow shock wave as the solution of the N-S Equations advances in time 
towards a steady state. 

Figure 5 shows the physical field, which constituted the 
computational domain and Figure 6 represents the transformed 5-n field 
used in the numerical computations. In this transformed field the body 
extends from C = 9 to ^ = 31 on the lower side. The upstream boundary 
transforms to the line on the top while the two vertical sides 
represent the downstream boundary. A fairly compact field with ^9 
lines in the ^-direction and 35 lines in the n-direction was used. 
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Even so the computer program required 62K of core capacity on ? DNIVAC 
1100 series computer. As the n-line spacing was already sparse to 
begin with, the scope for mesh refinement in the region of the shock 
was very much restricted. In Figure 7 the concentration of the rj-lines 
in the region of the shock is exhibited. 

Figures 8 and 9 are the Mach contour plots at the characteristic 
times of 0,8 and 3.2 and depict the progression of the bow shock wave 
from the body to its steady state stand-off distance. The Mach contour 
interval is 0.1. The sonic lines between the bow shock and the body, 
in which region the flow is wholly sub-sonic, are indicated in Figure 
9. It can be seen from this figure that aft of the body too, a sub- 
sonic region exists which extends up to a distance of about 2 times the 
diameter of the cylinder from its center. Behind this sub-sonic region, 
the flow again is super-sonic. In the field shotvn in Figure 5, the 
computational domain downstream of the body was limited by a semi- 
circle of radius 2.5 and the boundary conditions on this exit plane, 
as mentioned in Sect. 4.5, were established by solving the Euler's 
Equations on it. Since the do;<mstream boundary is located beyond the 
sub-sonic region and wholly in a super-sonic field of flow the use of 
the Euler's Equations is thus seen to be perfectly valid and accurate. 
In Figure 5 it can be noticed that the upstream boundary of the field 
starts and ends vertically above and below the clyinder respectively. 
This was dictated by the need of having to prescribe the free stream 
conditions at least up to those points. Figure 10 shows the Mach 
contours at time 3.2 for the field with mesh refinement. 

Figures 11 through 18 depict the variation of density, pressure, 
temperature and velocity from the front stagnation point of the cylinder 
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to the upstream boundary along the symmetry line (I = 20) in the steady 
state. In the physical field the front stagnation point is located at 
X = ~ 0.5 and the upstresua boundary* on the line of symmetry, at 
X = - 1.25, The absolute values of x are indicated in the plots. 

Figure 11 is the density distribution without mesh refinement, while 
12 is the density variation with mesh refinement in the shock region. 
While the trend of density, pressure, temperature and velocity distri“ 
butions seems satisfactory, the shock stand-off distance is more than 
what the ideal theories of Refs. [13] and [15] predict. This effect is 
due to the introduction of numerical dissipative terms, necessitated 
by the need to damp out wiggles. The disturbance of the smooth 
variation of the brought about by the mesh refinement seems 

to have given rise to the oscillations in the density, pressure, 
temperature and velocity profiles depicted in Figures 12, 14, 16, and 
18 respectively. Considering the coordinate system to be Independent 
of time could also have contributed to these oscillations. 

Figure 19 shows the variation of the coefficient of pressure (C ) 

P 

along the upper half of the cylinder from the front stagnation point to 
the trailing edge, while in Figure 20, is plotted for the field with 
mesh refinement. 

Figures 21 and 22 show the distribution of wall pressures from 0° 
to 180® normalized with the stagnation value without and with mesh 
refinement respectively. The mesh refinement, mild as it is, has not 
made any appreciable change in the pressure distribution on top of the 
cylinder as seen from Figures 20 and 22. In Figure 23 the results 
plotted in Figure 21 are compared with the experimental results quoted 
in Ref. [5] up to the 90“ point. It is seen that the numerical 



36 











solution of the compressible N~S Equations, for the free stream 
conditions considered, yields results which are quite in agreement with 
those obtained by experiments. 

Figures 24 and 25 show the computational domain and the pressure 
distribution on a cylinder with a fairing for = 4.6 and Re « 10. 

It is noticed that the pressure increases as we move closer to the 90® 
point, when compared with the results for the cylinder without a 
fairing. 




CHAPTER VI 


CONCLUSIONS 

6.1. In. general » In any numerical computational work, the finer the 
mesh system the more realistic the solution is expected to be, par- 
ticularly in regions where large gradients are encountered. This fact 
was amply made clear in the present research where it was found that 

as the intensity of the n-line spacing around the cylinder was increased 
the pressure distribution on it progressively approached the experi- 
mental results. 

6.2. It is also generally agreed upon that for any given intensity of 
the coordinate line spacing the farther the free stream boundary can be 
located from the body, the better the solution would be - subject to 
the limitations of computer time and storage. In the present investi- 
gation this factor was even more crucial as the bow shock tended to cut 
across the upstream boundary on which the free stream conditions were 
prescribed. Hence the upstream boundary was located as far upstream 
as we could and as parallel to the bow shock as possible, without 
compromising too much on the n-line spacing around the body. 

6.3. The influence of the dox^nstream flow conditions on the resultant 
flow field about a body in super-sonic flow seems to be a controversial 
topic. In this research for the same number of mesh points (39 x 35) 
various fields were considered in which the distance of the downstream 
boundary from the body varied, which in effect, varied the £-line 
spacing aft of the body. The best results were obtained for the 
boundary location that resulted in a finer ^-line spacing, thus 
conclusively proving that doxjustream flow conditions do affect the 
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numerical solution even in high Mach Number flows. This point is 
discussed in detail by Eoache [ 6 ]. 

6.4. The nearness of the downstream boundary to the body is not too 
detrimental so long as it is placed beyond the wake region aft of the 
body in which viscous effects do predominate and in which the flow is 
sub-sonic. The downstream boundary being located entirely in a super- 
sonic flow field, the Euler Equations may be used to establish accu- 
rately the downstream conditions. 

6.5. The wall density, except at the trailing edge, was calculated 
from the continuity equation written at the wall, using three-point 
forward difference approximations. However at the trailing edge, or 
the rear stagnation point, an explicit discretization of the continuity 
equation based on the leap-frog scheme proved beneficial in overcoming 
the negative densities that were otherwise occurring there. The same 
scheme when applied to the rest of the body, however, drove the wall 
densities down progressively, until negative values appeared in the 
neighborhood of the leading edge. 

6.6. Generally it is considered that a non-linear instability can 
result by using the average of the product (A.O.F.) scheme in 
approximating the convective derivatives, particularly in regions of 
flow reversal. In this research however, the A.O.F. scheme seemed to 
result in oscillations that were detrimental to the solution and thus 
the product of the average scheme was used. 

6.7. Normally we would expect a better solution by refining the mesh 
in the region of the shock where very high gradients of the flow vari- 
ables exist. In the present work mesh refinement in the shock region 










caused oscillations to build up quite rapidly there » while the pressure 
distribution on the body itself remained unaffected, l^hile no definite 
conclusions can be made at this stage, it is probable that the follow- 
ing factors caused the wiggles to appear in spite of the Shumann Filter. 

(a) The q-llne spacing in the computational domain was barely suf- 
ficient to start with; and the subsequent attraction of the q-lines 
into the region of the shock only caused a coarser mesh on either side 
of the shock, 

(b) Abrupt variations in the coordinate line spacing are not conductive 
to good results. It is possible that the mesh refinement that was 
attempted could have caused such uneven q-line spacing. 

(c) The solution could also be extremely sensitive to perturbations 
caused by the propagation of the shock by however small a distance, 
even though the mesh refinement was attempted at an almost steady-state 
stage of the solution. 

To sum up the work done, the numerical solution of the full N-S 
Equations for viscous compressible flows with a detached bow shock 
ahead of a two-dimensional circular cylinder has been obtained. The 
results agree quite well with experiment. Besides quite a steady state 
solution was obtained in a very short characteristic time. Before a 
conqiarison of the efficiency of the method can be made In relation to 
any other, the following points should be kept in mind. The compu- 
tational domain was small>^ 39 x 35 mesh points in all. Iterative 
convergence was obtained to a tolerance of for all variables at 

each time step (At = .01), in 6 iterations on an average. The total 
computer time required to achieve a steady state solution on an UMIVAC 
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1100 series computer, inclusive of the time required to generate the 
numerical curvilinear coordinate system, the coefficients of the 
Equations, etc,, is about 3 hours, A larger computational domain and 
a finer mesh will no doubt yield a better solution# The range of Mach 
Numbers and Reynolds Numbers for which this scheme is valid is yet to 
be tested. Nevertheless, judging from the results obtained, it Is 
believed that we are one more step ahead in using Computational Fluid 
Dynamics as a standard aerodynamic technique. 









i 




Up-stream 

boundary 


Down-stream 

boundary 


oi\’ic.d;al page is 

Oi?: PUOli QUALITY 


Figure 5. The Computational Domain 
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Figure 6* The Transformed Computational Domain 
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Figure 7. The Computational Domain with Mesh Refinement 
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Figure 9. Mach Contours at Time 3.2 
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Figure 11. Density Variation Across the Shock on Stagnation 
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Figure 19. Coefficient of Pressure Distribution on Top of 
Cylinder. 




























APPEHDIX A 


This Appendix susunarlzes the basic rules of tensor calculus » 

[16] used in transforming Eqs. (3.1) - (3.3), and the expanded form of 
the pertinent terms which appear in Eqs. (3.17) - (3.20). In all the 
formulae given below, repeated indices imply summation. 

Lex be the Cartesian coordinates and the general curvilinear 
coordinates. Then the metric coefficients are 


4 




95-" H 


J 


g 


ij _ kI ill 




ik .i 


g = det (gj^j) 


J “ 


3(Xj^,X2,X3) 

9(e^c^,c^) 


(A-1) 

(A-2) 

(A-3) 

(A-4) 

(A-5) 


The element of length ds is given by 

(ds)2 = 6^^ dx^ dXj 
= d?^ 


(A-6) 


Based on the above definitions, we collect other formulae from 
tensor calculus. 

Christoff el symbols: 



_ 1 11 
2 ® 


(■ 


9g£4 3g 


£k 


dx 


k 




!!ak) 


(A-7) 
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which Is symnietrlc in j»k. Contracting the i and j Indices, we have 




‘kr aC^ 


(A-! 


Covariant Derivative; 


,k 


. ^i r 
— T ■*■ r, V 
3 jk kr 


(A-l 


Divergence of a Vector; 


div V = v^) 


/E ac' 


(A-; 


Divergence of a Tensor jdelding Contravariant Components; 


, , . _»i 1 3 >. /- ikv , „i rs 

(div t) = — - ~-r C/g T ) + r T 

/g 31**^ 


(A- 


Laplacian of a Scalar 


35^ 3^ 


The transformation of Eqs. (3.1) - (3.3) is now direct, 
contravariant components, Eq. (3.1) becomes 


(A- 

Dsing 


|r + -p (P>^ v^) = 0 


(A- 


Equation (3.2) becomes 


gT (PV ) + — V 

^ 3r 


1 3 , ^ ik. , _i rs ^ 

(/g T )+r_T =0 


rs 


(A- 


where i = 1,2 for two dimensions 
Equation (3.3) becomes 

ST , 1 3 




(A- 


The expanded form of these equations are Eqs. (3.11) - (3.13). 
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The terms appearing in Eqs. (3.17) - (3.20) are 


e « i re iP _ - . p 


f OJET Oif P 

g - S 22 •** g <«22 + rf2> 


“ %2 ^^19 ^ 99 ) J + - (goo — - B 

12 12 22 ' ^ '•^22 95 8^2 9 ^^ 

8 D ^ 1 SD^ ^ 

' '“sT “ST + ■'u “" + 2ria + r^2°”> 


ap 


3 Px 


“ g ^^12 H ~ hi t ^8ii (ri2 + ^ 22 ^ 


- 812 O’}, + r.2 )} + £ (g 3K _ 3K. 




'J' = e {-^ (,/i vl) + ^ (^ v2) 


^“oW^^^%2f-B,2|^» 


Defining 


A iJL /„ 3 T ^ 

^®11 an ~ ®12 ag^ 


®ij ■ 

the terms ^ ^ 


V^ v 2 


Dll = 2 mI 6 ,, « - 


'22 «C - <5^2 “n ^®22 ^11 " G^2 ^12>‘^ 


(A-18) 


/A 1 rt \ 


* ^°22 ‘'12 - %2 '’ 22 >'’l 


(A-21) 


= u[(G,, u + vj - (u. + vj 


11 n 22 "12 '“C n" 


^®22 ^11 ■*■ ®11 ^12 “ *^12 ^^11 ^12^^“ 


’*' ^^22 ^12 ®11 ^22 " ®12 ^^12 ''' ^22^ 


22 - 


2p[G-« V — G_ „ v_ + (G« T « G- „ r?™)u 


11 n 12 "c '"11 ‘12 "12 *11" 


■*■^^11 ^22 " ®12 ^12^^^ 


K = x [u^ + + iTl^ + r^vl 


vl - Ku+ uCG^i G 22 u + G^2 ^22 V + u)vj^ 


+ h(G^ 2 ^22 “ ^‘^22^^^^^fl“ ^^®11 ®12 “ ■*■ ^®12^‘ 


- y{(G^2>^“ ®11 ®22 ^^2 


y2 - Kv - Gj ^2 


- m{(Gj^ 2>^“ ®12 ®22 ■ ''Ki 


+ w{(G^3^)2u + G^^ Gj^2 ^^^^2 


^ ’^^'^ll ^^12 “ + ^11 ®22 ^ "* ^^^^2 


(A-22) 

(A-23) 

•v-v}v^ 

» 

(A-24) 


(A-25] 



APPENDIX B 


EQUATION FOR NUFJERICAL COORDINATE GENERATION 

As in Appendix A, we shall denote the Cartesian coordinates as 

(the index i serving as a label only, having no tensor ial signif- 

i -i i 

icance) and the curvilinear coordiantes either as x , or x , or ? . 

It was mentioned in Section 2 that we need formulae in which the 
Cartesian coordinates x^ are treated as dependent variables, while the 
curvilinear coordinates C a-r® treated as independent variables. To 
achieve this goal, we shall use some formulae from general tensor 
calculus. In all the expressions given below, the repeated indices 
imply summation. 

Formula for the Second Derivative: 

Let X and x be two general coordinate systems. The formula for 
the s>-cond derivative [16] is 

9^x^ « t 3> ^ ft ^ 

3x^ 3?^ " 3^ ’ 3x^ 3?* (B-1) 


where P^. and r^, are the Christoff el symbols in the x^ and x^ coordl- 

jK jk 

nate systems respectively. (Refer to Eq. (A—7) for definition). Since 


we are considering the transformation between a Cartesian and a 
curvilinear system, hence either x^ or x^ is a Cartesian system. 

If x^ is a Cartesian system, then = 0. Writing x’" = x^ and 


x^ = 5^, we have 




3x 


pP 

35 *^ 3e“ *“ 35*’ 


(B-2) 


which is the formula for the second derivative of any Cartesian 
coordinate with respect to the curvilinear coordinates. 
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k 



Next, If x^ is a Cartesian system, then =0. Writing x^ = x. 

Am 1 

and = 5^, we have 


= _ r>r 3ii 


A> m St m 


(B-3) 


which is the formula for the second derivative of any curvilinear 
coordinate with respect to the Cartesian coordinates. 

The use of Eqs. (B-2) and (B-3) along with the equations 


„j.i . 9x. 

3€ _ gip £ 


3x„ 


3C' 


(B-4) 


a^P p 

yields a series of useful equations. 


(B-5) 


Inner nmltiplicatxon of Eq. (B-2) with and use of Eq. (B-5) 

9 X T ' 

r 


yields 


ij 


3g" ^'^s 

35^ 


(B-6) 


Using Eq. (B-4) in (B-7), we have 


3x 

rt s 


3^x 


ij 


35^ 


Introducing Eq, (B-4) in Eq, (B-3), we have 


= _ .r ip jq ^ ^ 


3x„ 3x 
£ m 


ij 


35^ 3^“ 

Another form of (B-8) can be obtained by using (B-7), which is 
_ _rt ^ip nq s 


rt ip jq 
3x. 3x = - 8 e ^ 


3x 3x„ 3x 
s £ m 

3C^ 3g^ 3g*^ 3g** 3g** 


(B-7) 


(B-8) 


(E-9) 


£ m 

It must be noted that the right hand sides of Eqs, (B-7) - (B-9) have 
differentiations with respect to as dssired. 
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Laplaclan: 

Setting A = in in Eq. (B-9) , using Eqs. (A-1) and (A-3) , 
over the index m, we get 




„2rr rt ij "s 
^ 5 ■ - g g — r 


S^x 


In two dimensions, writing 

- C, 5^ 


35^ 3£^ 


n, Xj^ - X, X2 


and using 


{which are a consequence of Eq. 

v2^ “ [(822 - 2g 


^12^®* g^^ ■ ®11^® 
(A“3)), we have 

12 »en *11 >'nn'*n 


- (8,2 *55 - 2gj^2 + 811 


3/2 


- t(822 852 - 28 i 2 + 8u %„>2{ 


- (822 ■ ^*12 * *11 


3/2 


where 


8 - gji 822 - <812)=' 


Similarly, using (B-3) and (B-4) , we easily obtain 

* “’u *n >’n • 2n * % >'e> ’’L *? 

where = C and = n. 

i 

The Laplaclan of a scalar function f(t ) is obtained as 


f). 


„2e 1 3 /-/"IKdr. 

- ~ — 7 C'^g g —r) 

✓g K H 


ik 3f 
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and summing 

(B-10) 

(B-11) 

(B-12) 


(B-13) 

(B-14) 

(B-15) 

(B-16) 

div (grad 
(B-17) 


rs 


which ill two dimensions has the form 


- HS22 " 283^2 ®11 ^Tin^ 


+ ( 2g^2 ^12 ” §22 ^11 “ ®11 ^ 22 ^ 


■*■ ^^§12 ^12 ®22 ^il “ §11 ^ 22 ^ 


(B-18) 


Unit Tangent and Normal Vectors: 

In the generation of coordinates we have taken clockwise traverse 
along the body contour as positive. Denoting the unit tangent and 
noannal vectors as t and n respectively, and k the constant unit vector 
normal to the plane of the curve, the vectors (t, n, k) are assumed to 
form a right-handed system. The unit tangent and normal vectors for 
the 5 - const, and n = const, cur'^es are 

1 




const . 


(ix + jy ) 

§22 


<5>n = const. “ 7^ JV 

§11 


(B-19) 


(B-20) 


(n)e _ (iy„ - jx„) 

C = const. ^ “-'ll n 


(n) 


§22 

1 


~ n = const. r- 
’^ll 


(- iyr + 4*r) 


(B-21) 

(B-22) 


where i and j are unit vectors along x and y respectively. The re- 
solved parts of the velocity vector along the cur/es 5 = const, and 
1 ) = const, are 


(v*t)_ 

- - C = const. 


j- ^“§12 '^§22^ 

•^8- 


(B-23) 


22 
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ff 










(B^26) 


where u and v are the contravariant con^jonents of the velocity vector 
V 9 which are related to the Cartesian components through Eq* (3*26). 
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APPENDIX C 


FINITE-DIFFERENCE APPROXIMATION OF THE N-S EQUATIONS 


This Appendix gives the details of the discretisation of the 
N-S Equations, Eq. (4.1), by a fully implicit difference approximation. 
The time derivatives are approximated by first-order backward differ- 
ences, while the spatial derivatives are approximated by central 
differences. The resulting system of finite-difference equations are 
solved by the S.O.R. method, - In what follows below the quantities at 
the previous time step are denoted by an overhead bar. Also 
and An do not occur in these equations as they both equal 1. 


Finite-differencing of equation (3.17): 


“ °i,J ^ 


■** ^®i,J+l ®i,J-l^^ 


In the S.O.R. scheme (C-1) is used as follows: 

Oi^j = (1-m) c^^jp + u,(o^^j)* 


(C-1) 


(C-2) 


In Equation (C-2) the superscript p denotes that the quantity is 
from the previous iteration. The second term on the R.H.S. of Equation 
(C-2) marked by an asterik is a provisional term defined by Equation 
(C-1), which la recovered from (C-2) for lo = 1. to = 0 indicates that 
no progress is being made in iterative convergence, which is a trivial 
case. 
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The finite-differencing of Equations (3-18) through (3.20) are 


carried out in a similar fashion; 


" ^1,J ■ 


+ ®i.J+l - ‘'l.J-l” 


- it [r^H + zri^N + rijM - el^^j 

A. ^ = (1-w) A, ,P + o)(A. T-)* 

Ijj XjJ 


" \.J - - ""i-l.j’ 




Atir^iH + 2F|2N + ^^i»3 


B = (i-(i) B. p + m(B y 

ijvJ X,J 




■*■ ^^i,J+l ^i,J-l^^ 


- At'F, . 
i,J 


(C- 

(C- 


(C- 

(C- 


h.J ■ «-“> 

The terms 0j 0 and Y that occur in Equations (C-3) , (C-5) and 


(C- 

(C- 


(C-7) are defined by Equations (A-16) , (A-17) and (A-18), repectively 
The finite-differencing of these equations is carried out next. 



Certain combinations of Christoffels and other terms that occur 
in groups frequently have been abbreviated as follows: 


DA - 1(^22 ' ^12^^12 




^ ^'^22^11 “ ‘^12^12^1,1 


DC “ 1^22^11 ~ ‘^12^12 ^i-l,j 
DD = [62/^2 “ Gl2^22Ji+l.j 


^®22^12 ~ ®12^22^i-l,j 


*-®22^il ®11^12 “ ®12^^il ■** ^12^^i,j 


- ^®22^il ■*■ ‘^11^12 ‘^12^^11 ^12^^i,j-l 


!/■ 

' : ' 


^^22^12 ^ll^L ~ ^12^^12 ^22^^i,j+l 


^®22^12 ’*■ ®11^22 “ ®12^^12 ^22^^i,j-l 


^^^11^2 - ^12^il^i,j 


^^22^12 ■ ®12^22^i,j 


*•^22^12 ^ll^L " ®12^^i? ‘*' ^22^^i,j 


*•^1^22 " ^12^12^ i.j 

*•^22^11 ^11^12 “ ‘^12^^11 '‘' ^12^^i-l,j 


^^2^12 ■*■ ^11^22 " ^2^^i2‘^^22^U+l,j 
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tS2^12 - ^12^"i2 "’225Ji-l,J 



^®11^12 " ®12^il^i,j+l 
^®11^12 “ ®12^il^i,j-l 


^®11^22 “ *^12^12^ i,j+l 


^^11^22 ■" ®12^12^i,j-l 


With these abbreviations the finite-difference approximation of 
Equation (A-16) takes the following shape. 

®i,j " ^‘^12^^’^^i,j ^Pi,j+1 " Pi.j-1^ 

+ £( 622 / 2 )^^^ - “ij) 


^^i,j ^i-l,j^^'^i,j “ "i-l,j^^ 


+ £( 622 / 4 ) ^"’'i+lj+l ^i+l,j-l^ 


" ^i-l,j^Vl,j+l " VlJ-l^^ 


+ XG 22 / 2 ) [((rj,fr|2) Xu),^,^. 


- ((rli + r|2) Xu)^_^ 


+ £(622/2) E((r^2 ■*■ ^22^^^^i+l,j 


(Cri2 + ry xv)^_ 


^^^12^^^^ ^^i,j+l ^“i+l,j+l " '^i-ljj+l^ 


\,j-l ^“i+l,j-l " Vl,j-1^^ 


'*' ^i, j+1^ j+1 " "^i.j^ 


e(Gj^2/2) + f22) Xu)^ - ((rj, + T?J Xu) 
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- b(G^ 2/2) + r|2)Xv), - ((ri2 

+ e[((pG22)^^j + ' “i,j^ 


((U 622 )ij + 


(e/2)[(ii6j^)^^^^j (“i+i,j+l " "i+l,j-l^ 




+ e[DB(pu)^^^^^ - DC(pu)^^^^^] 


+ e[DD(yv)^^3_^^ - DE(uv)^_3_^ . 3 


+ (e/2)[((u6j^^)^ . + CuG3^i)i,j+lH“i,j+l “ «i,j) 

+ Ce/4) E(uG22)i^j+l - Vl,j+1^ 




- (e/4) tuGj^ 2 >i,j+l ^"i+l,j+l " Vl.j+1^ 


~ ^'^®12^k,j-l ^“i+l,j-l “ “i-l,j-l 


)3 


(c/2) [((pGj^ 2 ^i,j ■*■ 




+ (e/2) [DF(uu).^^^^ - DG(yu)^^^_^] 


+ (e/2) [DH(pv)^^^^^ - DICmv)^^._^1 
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i-J J J J 1 


L-. 


_1 :.l 


+ BB + DF + dj) 




K • -1 • T *• tlj -1 • l)l 

i,J-l 1+1, 3-1 i-1,3-1'’* 


+ (e/2)(G -). [(A. . + A ..,)(v. , , - v. .) 

11 1,3 1,3 1,J+1 1,1+1 lij 


- (^^ ^ + ^ i)(v,. . “V . ,) 


- (e/2)(6j^2^i,3 ^i+l,j^ ^"i+l,j “ “l.j^ 


“ ^^i ,3 ■*■ ■ Vl,j^^ 


(E/4)(G^2>i,j ^^i+l,j ^^i+l,j+l ~ ’^i+1,3-1^ 


! 


■*■ ^'"i, 3 ^^i+l,j " ’^i- 1 , 3 ^^® 


^^±,j ^^ 1 , 3+1 " ’^ 1 , 3 - 1 ^ *-^ 12*^12 “ ^L®ll^i,j 


m- ■ . -; 

T' 

II'- i 


+ 2e(pv)^^^ dk + dl + dm] 

‘‘'ij “ <® 12 ^’^>i,j ^Pi+ 1,3 “ 

" ^®ll^‘^^i,j ^Pi,j +1 " ^i, 3 -l^ 


(C-9) 


+ ((rj^ + + i’f 2 ^^“^i,j-l 


+ (( r ]_2 + ^ 22 ^ ^^^ 1 , 3+1 " ^^’^12 ■*' ^ 22 ^^^^ i , j - l ^ 


^i-1,3 ^^i-1,3+1 ” ^1-1, j-1^^ 
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- (e/2) (63^)^^ j **■ 

^(^11 **■ ■*’ ^22^^^^i+l,j 


- (Cri2 + q2>^^h-i,2 

+ (e/4) ^“i+l,j+l “ ^i+l,j-l^ 




+ (e/2) [((V‘G 22 )i^j + ^^i+l,j “ ^i,j^ 


((jiGi2)^^j + (l^‘^i2^i+l,j^ ^“l+l,j ” 




- (e/4) [(yG3^2^i+l,j ^"^i+ljj+l ~ '^i+ljj-l^ 


- (yG^2)i-i,j ("^i-lj+l “ 


^ (e/2) [(yu).^^^. DF - (yu)^^,^. DO 




+ e[((y6j^^)^_^ + (^Gii)ij+l5(^l,j+l ~ ^i,j- 
((yGj^j^)i^^ + (^‘^ll^i,j-l^^^i,j " 

- (e/2)[(yGj^2'i,j+i ^^i+l,j+l “ ^i-l,j+l^ 
(yGj^2^i,j-l (\+l,j-l “ 
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'It* 










+ (liv).^ BT - (PV)^^ DD] 

■*' ^^i,j^^”i+l,j “ “i-li,j^^^ilS2^i,j " ^^12^12^1, 

■*■ " '^i-l,j^^^^l2®22^i,j " ^^22®12^i,j^ 

■*■ ^‘^i,i+l “ '^i,j-1^^^^22®ll^i,j ” ^^12^12^1,3^^ 

+ 2e(,u)^^^[(rf^)^^. DB + (rj^).^. DF t- iVy^^. DJ] 

+ 2e(yv)^^. [(Vy^^. DK + . DL + DM] 

+ (a^/2)[((uG22)j^^j + ^’^®22^i+l,j^ j " ^i,j^ 


^^^® 22 ^i,j ■*' ^^® 22 ^i-l,j^^'^i,j “ 


(a^/4)[(M6^2^i+lj ^’^i+i,j+i " \+l,j-l^ 


(pG^ 2 ^i-l,j 


+ (a^/2)l((iiG^^).^^ + - T. .) 




(a^/4)[(pG^2^ij+i *^'’^1+1, j+l “ "*^1-1, j+1^ 


^*'® 12 ^i,j-l ^^i+l,j-l “ ^ 1 - 1 , j- 1 ^^ 


(C-10) 


(C-11) 
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AFTENDIX-D 


The non-diraensionalization of the N-S Equations is given in many 
standard text books. The non-dimensionalization of the other equations 
such as the Sutherland viscosity law, the thermod 3 mamic equations of 
state, etc., are given in this Appendix. 

Even in the range of super-sonic Mach Numbers, M < 6 in air, the 
temperature rise in the gaseous stream is high enough to force us to 
take into account the effect of temperature on the properties of the 
gas, in particular, on its viscosity. The kinematic viscosity of most 
gases, and of air among them, increases considerably as the temperature 
is increased. The temperature dependence of viscosity p(T) must be 
obtained from experiments. The Sutherland viscosity law is; 

V* = C[(T*)^-^/{T* + S*)] (D-1) 

where C and are constants and for air have the values C = 2.27 x 10“® 
and S* = 110“K. Dividing Equation (D-1) by the corresponding ex- 
pression for free stream conditions, we get; 

(y*/p*) = (T*/T*)^*^ I(T* + S*)/(T* + S*)] 

to CO to 

T 

or in nondimensional terms 

P = T^-^ (1 + S^)/(T + S^) (D-2) 

where S, = S*/T* 

1 » 

The nondimensional equation of state is obtained as follows; 
sxnce p = p (c - c )T 

^ p V 


(D-3) 


t 




} 

i 




. »(c^ - v(i*/»;^)(Cp/=j,) 

= pCy - 1) C T*T/(y U*2) 

= p(Y^ - 1) h*T/(Y^U*2) 

= pCYj, - 1) T/(Yp(Yj, - 1)M^) 

i.e. f the nondlmenslonal equation of state is: 

P = PT/Y^M2 

Nondimenslonalization of the equation for temperature, now 
4' = p*e* + (l/2)p*lv*P , thus 

y = 'i'*/p*U*2 + p[(e*/U*^) + lyl^/2], or 

('i'/p) = I(c^T*/U*^) + jvlV2] 

“ [(%TT*Cp/CpD*2) + lvP/2] 

= [(TT%p/Y^U*2) + lvi2/2] 

= UT/y^CYo - 1)m 2) + lvP/2, or 

T “ DM^I'I'/P) - IVp/Z] 

o o ^ 

where lyj^ = g vV 

and are the contravariant components, u and v, of V. Thus 
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(D-4) 


(D-5) 


(D-6) 


The local Mach Number , which was used in plotting the Mach 


contours, is obtained as follows: 

= lYl/a = 1 v1//(y^p/o) , or 


- 7 - 


(2 
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